###################################################
# R SCRIPT
# 
#             by  : Kenneth Bunker (Universidad San Sebastián)
#                 : Gabriel Negretto (Universidad Carlos III de Madrid)
# corresping      : kenneth.bunker@uss.cl
# journal         : Party Politics
# article         : The Party System Effects of 
#                   Unstable Electoral Rules in Latin America
#
# CITE DOI    : Bunker, Kenneth; Negretto, Gabriel, 2023, "Replication 
#                 Data for: The Party System Effects of Unstable Electoral 
#                 Rules in Latin America", https://doi.org/10.7910/DVN/QPA1T3, 
#                 Harvard Dataverse  
#
#   version: RStudio 2023.06.1+524 "Mountain Hydrangea" Release 
#   (547dcf861cac0253a8abb52c135e44e02ba407a1, 2023-07-06) for macOS
#   Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_7) AppleWebKit/537.36 
#   (KHTML, like Gecko) RStudio/2023.06.1+524 Chrome/110.0.5481.208 
#   Electron/23.3.0 Safari/537.36
#
###################################################

## Drop everything and clean
rm(list=ls())

###################################################
# WORKING DIRECTORY
###################################################

setwd("/Users/kennethbunker/Dropbox/1. Negretto - Revisiting interactions/2023")

###################################################
# OPEN DATASETS
###################################################

library(dplyr)
library(robustbase) 
library(foreign)
library(sandwich)
library(lmtest)
library(estimatr)
library(readxl)
library(margins)
library(ggplot2)
library(interactions)
library(jtools)
library("interplot")
library("estimatr")
library(vtable)
library(huxtable)
library(stargazer)
library(rootSolve)

###################################################
# OPEN DATASET
#
#
###################################################

bunkernegretto2022 <- read_excel("../Data/replication data (HD).xlsx")

##############################################
# Choose measurements
#
#
###################################################

## ethnic
bunkernegretto2022$ethnic <- bunkernegretto2022$e2
# proximity
bunkernegretto2022$proximity <- bunkernegretto2022$p1

###################################################
# Generate variables
#
#
###################################################

## main interaction
bunkernegretto2022$logdm <- log(bunkernegretto2022$dep_dm_1tier)
## other interaction
bunkernegretto2022$sp <- bunkernegretto2022$dep_dm*bunkernegretto2022$dep_as
## other
bunkernegretto2022$legacy <- bunkernegretto2022$number+bunkernegretto2022$spell_years_elecc

###################################################
# TABLE 1: DESCRIPTIVE STATS
#
#
###################################################

## descriptive by country
st(bunkernegretto2022, vars = c("enpv_bn", "ethnic", "logdm", 
                                "dep_dm_1tier", "dep_as_1tier",
                                "dep_elec_regime", "dep_elec_acc",
                                "upper", "enc_bn", "proximity", 
                                "time", "number", "stability"), group="country_name", group.long = TRUE)

###################################################
# TABLE 2/MODEL 1: INTERACTIVE - ALL ELECTIONS
#
#
###################################################

m1 <- bunkernegretto2022
model1 <- lm(enpv_bn ~ ethnic + logdm + (ethnic*logdm) + 
               upper + enc_bn + proximity + 
               (ethnic*upper) + (enc_bn*proximity) + 
               time, data=m1)
# Stargazer output (with and without RSE)
cl_robust1 <- coeftest(model1, vcov = vcovCL, type = "HC1", cluster = ~ country_name)
se_robust1 <- cl_robust1[, 2]
cl_robust1

###################################################
# TABLE 2/MODEL 2: INTERACTIVE - EARLY ELECTIONS
#
#
###################################################

m2 <- subset(bunkernegretto2022, number<6)
model2 <- lm(enpv_bn ~ ethnic + logdm + 
               (ethnic*logdm) + 
               upper + enc_bn + proximity + 
               (ethnic*upper) + (enc_bn*proximity) + 
               time, data=m2)
# Stargazer output (with and without RSE)
cl_robust2 <- coeftest(model2, vcov = vcovCL, type = "HC1", cluster = ~ country_name)
se_robust2 <- cl_robust2[, 2]
cl_robust2

###################################################
# TABLE 2/MODEL 3: INTERACTIVE - LATE ELECTIONS
#
#
###################################################

m3 <- subset(bunkernegretto2022, number>5)
model3 <- lm(enpv_bn ~ ethnic + logdm + 
               (ethnic*logdm) + 
               upper + enc_bn + proximity + 
               (ethnic*upper) + (enc_bn*proximity) + 
               time, data=m3)
# Stargazer output (with and without RSE)
cl_robust3 <- coeftest(model3, vcov = vcovCL, type = "HC1", cluster = ~ country_name)
se_robust3 <- cl_robust3[, 2]
cl_robust3

###################################################
# FIGURE 1
#
#
###################################################

## TABLE 1: MODEL 3
m3 <- subset(bunkernegretto2022, number>5)
figure1 <- lm_robust(enpv_bn ~ ethnic + logdm + (ethnic*logdm) +
                       upper + enc_bn + proximity +
                       (ethnic*upper) + (enc_bn*proximity) +
                       time, data=m3,
                     clusters=country_name,
                     se_type="stata")
summary(figure1)
## x-axis
vertical <- seq(0,6,.1)
# get betas
b0 <- figure1$coefficients[1] 
b1 <- figure1$coefficients[2]  
b2 <- figure1$coefficients[3] 
b3 <- figure1$coefficients[8]  
# get variance
v0 <- vcov(figure1)[1,1] 
v1 <- vcov(figure1)[2,2]  
v2 <- vcov(figure1)[3,3]  
v3 <- vcov(figure1)[8,8]  
## covariance
covb1 <- vcov(figure1)[2,8]  
covb2 <- vcov(figure1)[3,8] 

## Figure 1
pdf(file=paste("../Party Politics/graphics/Figure 1.pdf", sep=""), width = 6, height = 6)
z1=0  #supply lower bound for z
z2=5  #supply upper bound for z
ci <- 1.684 #90%
z <- seq(z1,z2,.1)
fz <- c(z,z)
fz <- c(z,z)
y1 <- (b1+b3*z)+(ci*sqrt(v1+(2*z*covb1)+((z^2)*v3)))
y2 <- (b1+b3*z)-(ci*sqrt(v1+(2*z*covb1)+((z^2)*v3)))
fy <- c(y1,y2)
fline <- (b1+b3*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,
     xlab='Moderator: DM (Log)',
     ylab='Simple Slope: Ethnic Heterogeneity',
     main="Figure 1. The marginal effect of Ethnic Heterogeneity on \nthe Effective Number of Electoral Parties \nby Average District Magnitude",
     xlim=c(0,5))
lines(z,fline)
f0 <- array(0,c(length(z)))
## other
lines(z,f0,col=8)
zerocross <- round(uniroot.all(approxfun(z,y2), interval = range(c(0,10))),digits=3) ## check for zero-crossings 
abline(v=zerocross,col=4,lty=2)
legend("topright", cex=.5,  bty = "n", bg="white", legend = paste0("Significant above ",zerocross))
legend("bottomright", cex=.5,  bty = "n", bg="white", legend = "Table 2, Model 3")
dev.off()

###################################################
# TABLE 3
# 
# 
###################################################

table3rb <- lm_robust(enpv_bn ~ ethnic + logdm +  
                        number +
                        (logdm*ethnic) +
                        (logdm*number) +
                        (number *ethnic) +
                        (logdm*number*ethnic) +
                        upper + 
                        (upper*ethnic) + 
                        proximity + 
                        enc_bn + 
                        (proximity*enc_bn) + 
                        time,
                      data=bunkernegretto2022,
                      clusters=country_name,
                      se_type="stata")
summary(table3rb)

###################################################
# FIGURE 2: TWO (2) ELECTIONS
# 
#
###################################################

table3rb <- lm_robust(enpv_bn ~ ethnic + logdm +  
                        number +
                        (logdm*ethnic) +
                        (logdm*number) +
                        (number *ethnic) +
                        #(number*ethnic) +
                        (logdm*number*ethnic) +
                        upper + 
                        (upper*ethnic) + 
                        proximity + 
                        enc_bn + 
                        (proximity*enc_bn) + 
                        time,
                      data=bunkernegretto2022,
                      clusters=country_name,
                      se_type="stata")
summary(table3rb)

## EXTRACT VALUES
## x-axis
vertical <- seq(0,6,.1)
# get betas
b0 <- table3rb$coefficients[1]  
b1 <- table3rb$coefficients[2]  
b2 <- table3rb$coefficients[3]  
b3 <- table3rb$coefficients[4]  
b4 <- table3rb$coefficients[9]  
b5 <- table3rb$coefficients[10] 
b6 <- table3rb$coefficients[11] 
b7 <- table3rb$coefficients[14] 
# get variance
v0 <- vcov(table3rb)[1,1]    
v1 <- vcov(table3rb)[2,2]    
v2 <- vcov(table3rb)[3,3]    
v3 <- vcov(table3rb)[4,4]    
v4 <- vcov(table3rb)[9,9]    
v5 <- vcov(table3rb)[10,10]  
v6 <- vcov(table3rb)[11,11]  
v7  <- vcov(table3rb)[14,14] 
## covariance
covb4_b1 <- vcov(table3rb)[9,2]   
covb5_b1 <- vcov(table3rb)[10,2]  
covb7_b1 <- vcov(table3rb)[14,2]  
covb5_b4 <- vcov(table3rb)[10,9]  
covb7_b4 <- vcov(table3rb)[14,9]  
covb7_b5 <- vcov(table3rb)[14,10] 
covb2_b0 <- vcov(table3rb)[3,1]   
covb3_b0 <- vcov(table3rb)[4,1]   
covb6_b0 <- vcov(table3rb)[11,1]  
covb3_b2 <- vcov(table3rb)[4,3]   
covb6_b2 <- vcov(table3rb)[11,3]  
covb6_b3 <- vcov(table3rb)[11,4]

## number=2
pdf(file=paste("../Party Politics/graphics/Figure 2.pdf", sep=""), width = 6, height = 6)
z1=0  #supply lower bound for z
z2=5   #supply upper bound for z
ci <- 1.684 #90%
nelections <- 2
z <- seq(z1,6,.1)
fz <- c(z,z)
y1 <- (b1+b4*z+b5*nelections+b7*z*nelections)+ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
y2 <- (b1+b4*z+b5*nelections+b7*z*nelections)-ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
fy <- c(y1,y2)
fline <- (b1+b4*z+b5*nelections+b7*z*nelections)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,
     xlab='Moderator: District Magnitude (Log)',
     ylab='Simple Slope: Ethnic Heterogeneity',
     main='Figure 2. Marginal Effect of Ethnic Heterogeneity on ENPV, \nby Average District Magnitude when Number=2',
     xlim=c(0,5))
lines(z,fline)
f0 <- array(0,c(length(z)))
lines(z,f0,col=8)
zerocross <- round(uniroot.all(approxfun(z,y2), interval = range(c(0,10))),digits=3) ## check for zero-crossings 
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)
legend("bottomright", cex=.5,  bty = "n", bg="white", legend = "Table 3")
dev.off()

###################################################
# FIGURE 3: NINE (9) ELECTIONS
# 
#
###################################################

table3rb <- lm_robust(enpv_bn ~ ethnic + logdm +  
                        number +
                        (logdm*ethnic) +
                        (logdm*number) +
                        (number *ethnic) +
                        (number*ethnic) +
                        (logdm*number*ethnic) +
                        upper + 
                        (upper*ethnic) + 
                        proximity + 
                        enc_bn + 
                        (proximity*enc_bn) + 
                        time,
                      data=bunkernegretto2022,
                      clusters=country_name,
                      se_type="stata")
summary(table3rb)

# get betas
b0 <- table3rb$coefficients[1]
b1 <- table3rb$coefficients[2] 
b2 <- table3rb$coefficients[3]
b3 <- table3rb$coefficients[4]
b4 <- table3rb$coefficients[9] 
b5 <- table3rb$coefficients[10]
b6 <- table3rb$coefficients[11]
b7 <- table3rb$coefficients[14]
# get variance
v0 <- vcov(table3rb)[1,1]
v1 <- vcov(table3rb)[2,2]
v2 <- vcov(table3rb)[3,3]
v3 <- vcov(table3rb)[4,4]
v4 <- vcov(table3rb)[9,9] 
v5 <- vcov(table3rb)[10,10]
v6 <- vcov(table3rb)[11,11]
v7  <- vcov(table3rb)[14,14]
## covariance
covb4_b1 <- vcov(table3rb)[9,2]
covb5_b1 <- vcov(table3rb)[10,2]
covb7_b1 <- vcov(table3rb)[14,2]
covb5_b4 <- vcov(table3rb)[10,9] 
covb7_b4 <- vcov(table3rb)[14,9]
covb7_b5 <- vcov(table3rb)[14,10]
covb2_b0 <- vcov(table3rb)[3,1]
covb3_b0 <- vcov(table3rb)[4,1]
covb6_b0 <- vcov(table3rb)[11,1]
covb3_b2 <- vcov(table3rb)[4,3] 
covb6_b2 <- vcov(table3rb)[11,3]
covb6_b3 <- vcov(table3rb)[11,4]

## number=9
pdf(file=paste("../Party Politics/graphics/Figure 3.pdf", sep=""), width = 6, height = 6)
z1=0  #supply lower bound for z
z2=5   #supply upper bound for z
ci <- 1.684 #90%
nelections <- 9
z <- seq(z1,6,.1)
fz <- c(z,z)
y1 <- (b1+b4*z+b5*nelections+b7*z*nelections)+ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
y2 <- (b1+b4*z+b5*nelections+b7*z*nelections)-ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
fy <- c(y1,y2)
fline <- (b1+b4*z+b5*nelections+b7*z*nelections)
plot(fz,fy,type='p',pch=".", cex=1, font=2,font.lab=2,col=2,
     xlab='Moderator: District Magnitude (Log)',
     ylab='Simple Slope: Ethnic Heterogeneity',
     main='Figure 3. Marginal Effect of Ethnic Heterogeneity on ENPV, \nby Average District when Number=9',
     xlim=c(0,5))
lines(z,fline)
f0 <- array(0,c(length(z)))
lines(z,f0,col=8)
zerocross <- round(uniroot.all(approxfun(z,y2), interval = range(c(0,10))),digits=3) ## check for zero-crossings 
abline(v=zerocross[[1]],col=4,lty=2)
abline(v=zerocross[[2]],col=4,lty=2)
legend("topright", cex=.5,  bty = "n", bg="white", legend = paste0("Significant between ",zerocross[[1]]," and ",zerocross[[2]]))
legend("bottomright", cex=.5,  bty = "n", bg="white", legend = "Table 3")
dev.off()

###################################################
# FIGURE A1: ALTERNATIVE PLUG-IN ELECTIONS
# 
#
###################################################

# ## number=variable (plug in): 
# nelections <- 12 ##        <------------------- plug-in number of elections
# pdf(file=paste("../Party Politics/graphics/Figure A1 (",nelections," elections).pdf", sep=""), width = 6, height = 6)
# z1=0  #supply lower bound for z
# z2=5   #supply upper bound for z
# ci <- 1.684 #90%
# z <- seq(z1,6,.1)
# fz <- c(z,z)
# y1 <- (b1+b4*z+b5*nelections+b7*z*nelections)+ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
# y2 <- (b1+b4*z+b5*nelections+b7*z*nelections)-ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
# fy <- c(y1,y2)
# fline <- (b1+b4*z+b5*nelections+b7*z*nelections)
# plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,
#      xlab='Moderator: District Magnitude (Log)',
#      ylab='Simple Slope: Ethnic Heterogeneity',
#      main=paste0("Number=",nelections),
#      xlim=c(0,5))
# lines(z,fline)
# f0 <- array(0,c(length(z)))
# lines(z,f0,col=8)
# legend("bottomright", cex=.5,  bty = "n", bg="white", legend = "Table 3")
# dev.off()

###################################################
# TABLE 4, MODEL 1
# 
# 
###################################################

d_table4m1 <- subset(bunkernegretto2022, dep_elec_acc>2)
table4m1 <- lm_robust(enpv_bn ~ ethnic + logdm +  
                        (logdm*ethnic) +
                        upper + 
                        (upper*ethnic) + 
                        proximity + 
                        enc_bn + 
                        (proximity*enc_bn) + 
                        time,
                      data=d_table4m1,
                      clusters=country_name,
                      se_type="stata")
summary(table4m1)

###################################################
# TABLE 4, MODEL 2
# 
# 
###################################################

d_table4m2 <- subset(bunkernegretto2022, number < 6 & dep_elec_acc>2)
table4m2 <- lm_robust(enpv_bn ~ ethnic + logdm +  
                        (logdm*ethnic) +
                        upper + 
                        (upper*ethnic) + 
                        proximity + 
                        enc_bn + 
                        (proximity*enc_bn) + 
                        time,
                      data=d_table4m2,
                      clusters=country_name,
                      se_type="stata")
summary(table4m2)

###################################################
# TABLE 4, MODEL 3
# 
# 
###################################################

d_table4m3 <- subset(bunkernegretto2022, number > 5 & dep_elec_acc>2)
table4m3 <- lm_robust(enpv_bn ~ ethnic + logdm +  
                        (logdm*ethnic) +
                        upper + 
                        (upper*ethnic) + 
                        proximity + 
                        enc_bn + 
                        (proximity*enc_bn) + 
                        time,
                      data=d_table4m3,
                      clusters=country_name,
                      se_type="stata")
summary(table4m3)

###################################################
# FIGURE 4: TABLE 4, MODEL 3
# 
#
###################################################

d_table4 <- subset(bunkernegretto2022, number > 5 & dep_elec_acc>2)
figure4 <- lm_robust(enpv_bn ~ ethnic + logdm +  
                       (logdm*ethnic) +
                       upper + 
                       (upper*ethnic) + 
                       proximity + 
                       enc_bn + 
                       (proximity*enc_bn) + 
                       time,
                     data=d_table4,
                     clusters=country_name,
                     se_type="stata")
summary(figure4)
# get betas
b0 <- figure4$coefficients[1]
b1 <- figure4$coefficients[2]
b2 <- figure4$coefficients[3] 
b3 <- figure4$coefficients[8] 
# get variance
v0 <- vcov(figure4)[1,1]  
v1 <- vcov(figure4)[2,2] 
v2 <- vcov(figure4)[3,3]  
v3 <- vcov(figure4)[8,8]  
## covariance
covb1 <- vcov(figure4)[2,8]  
covb2 <- vcov(figure4)[3,8]  

## new
pdf(file=paste("../Party Politics/graphics/Figure 4.pdf", sep=""), width = 6, height = 6)
z1=0  #supply lower bound for z
z2=6  #supply upper bound for z
ci <- 1.684 #instead of 1.8331
z <- seq(z1,6,.1)
fz <- c(z,z)
fz <- c(z,z)
y1 <- (b1+b3*z)+(ci*sqrt(v1+(2*z*covb1)+((z^2)*v3)))
y2 <- (b1+b3*z)-(ci*sqrt(v1+(2*z*covb1)+((z^2)*v3)))
fy <- c(y1,y2)
fline <- (b1+b3*z)
plot(fz,fy,type='p', pch='.',font=2,font.lab=2,col=2,
     xlab='Moderator: DM (Log)',
     ylab='Simple Slope: Ethnic Heterogeneity',
     main="Figure 4. The marginal effect of Ethnic Heterogeneity on \nthe Effective Number of Electoral Parties \nby Average District Magnitude",
     xlim=c(0,5))
lines(z,fline)
f0 <- array(0,c(length(z)))
lines(vertical,f0,col=8)
zerocross <- round(uniroot.all(approxfun(z,y2), interval = range(c(0,10))),digits=3) ## check for zero-crossings 
abline(v=zerocross,col=4,lty=2)
legend("topright", cex=.5,  bty = "n", bg="white", legend = paste0("Significant above ",zerocross))
legend("bottomright", cex=.5,  bty = "n", bg="white", legend = "Table 4, Model 3")
dev.off()

###################################################
# FIGURE 5: TABLE 4, MODEL 3
# 
# 
###################################################

d_table4 <- subset(bunkernegretto2022, number > 5 & dep_elec_acc>2)
figure5 <- lm_robust(enpv_bn ~ proximity + 
                       enc_bn + 
                       (proximity*enc_bn) + 
                       ethnic + logdm +  
                       (logdm*ethnic) +
                       upper + 
                       (upper*ethnic) + 
                       proximity + 
                       enc_bn + 
                       time,
                     data=d_table4,
                     clusters=country_name,
                     se_type="stata")
summary(figure5)
# get betas
b0 <- figure5$coefficients[1]
b1 <- figure5$coefficients[2]
b2 <- figure5$coefficients[3] 
b3 <- figure5$coefficients[8] 
# get variance
v0 <- vcov(figure5)[1,1]  
v1 <- vcov(figure5)[2,2] 
v2 <- vcov(figure5)[3,3]  
v3 <- vcov(figure5)[8,8]  
## covariance
covb1 <- vcov(figure5)[2,8]  
covb2 <- vcov(figure5)[3,8]  

## new
pdf(file=paste("../Party Politics/graphics/Figure 5.pdf", sep=""), width = 6, height = 6)
z1=0  #supply lower bound for z
z2=10  #supply upper bound for z
ci <- 1.684 #instead of 1.8331
z <- seq(z1,10,.1)
fz <- c(z,z)
y1 <- (b1+b3*z)+(ci*sqrt(v1+(2*z*covb1)+((z^2)*v3)))
y2 <- (b1+b3*z)-(ci*sqrt(v1+(2*z*covb1)+((z^2)*v3)))
fy <- c(y1,y2)
fline <- (b1+b3*z)
plot(fz,fy,type='p', pch='.',font=2,font.lab=2,col=2,
     xlab='Moderator: Effective Number of Presidential Candidates',
     ylab='Simple Slope: Proximity',
     main="Figure 5. The marginal effect of Proximity on \nthe Effective Number of Electoral Parties \nby Number of Presidential Candidates",
     xlim=c(0,10))
lines(z,fline)
f0 <- array(0,c(length(z)))
lines(z,f0,col=8)
zerocross <- round(uniroot.all(approxfun(z,y1), interval = range(c(0,10))),digits=3) ## check for zero-crossings 
abline(v=zerocross,col=4,lty=2)
legend("topright", cex=.5,  bty = "n", bg="white", legend = paste0("Significant below ",zerocross))
legend("bottomright", cex=.5,  bty = "n", bg="white", legend = "Table 4, Model 3")
dev.off()

###################################################
# TABLE 5
# 
#
###################################################

table5rb <- lm_robust(enpv_bn ~ ethnic + logdm +  
                        dep_elec_acc +
                        (logdm*ethnic) +
                        (logdm*dep_elec_acc) +
                        (dep_elec_acc *ethnic) +
                        (dep_elec_acc*ethnic) +
                        (logdm*dep_elec_acc*ethnic) +
                        upper + 
                        (upper*ethnic) + 
                        proximity + 
                        enc_bn + 
                        (proximity*enc_bn) + 
                        time,
                      data=bunkernegretto2022,
                      clusters=country_name,
                      se_type="stata")
summary(table5rb)

###################################################
# FIGURE 6: TWO (2) ELECTIONS
# 
#
###################################################

table5rb <- lm_robust(enpv_bn ~ ethnic + logdm +  
                        dep_elec_acc +
                        (logdm*ethnic) +
                        (logdm*dep_elec_acc) +
                        (dep_elec_acc *ethnic) +
                        (dep_elec_acc*ethnic) +
                        (logdm*dep_elec_acc*ethnic) +
                        upper + 
                        (upper*ethnic) + 
                        proximity + 
                        enc_bn + 
                        (proximity*enc_bn) + 
                        time,
                      data=bunkernegretto2022,
                      clusters=country_name,
                      se_type="stata")
summary(table5rb)

# get betas
b0  <- table5rb$coefficients[1]   # (intercept)
b1  <- table5rb$coefficients[2]   # ethnic
b2  <- table5rb$coefficients[3]   # logdm
b3  <- table5rb$coefficients[4]   # dep_elec_acc
b4  <- table5rb$coefficients[9]   # (ethnic:logdm)
b5  <- table5rb$coefficients[10]  # (logdm:dep_elec_acc)
b6  <- table5rb$coefficients[11]  # (ethnic:dep_elec_acc)
b7  <- table5rb$coefficients[14]  # (ethnic:logdm:dep_elec_acc)
# get variance
v0 <- vcov(table5rb)[1,1]
v1 <- vcov(table5rb)[2,2]
v2 <- vcov(table5rb)[3,3]
v3 <- vcov(table5rb)[4,4]
v4 <- vcov(table5rb)[9,9] 
v5 <- vcov(table5rb)[10,10]
v6 <- vcov(table5rb)[11,11]
v7  <- vcov(table5rb)[14,14]
## covariance
covb4_b1 <- vcov(table5rb)[9,2]
covb5_b1 <- vcov(table5rb)[10,2]
covb7_b1 <- vcov(table5rb)[14,2]
covb5_b4 <- vcov(table5rb)[10,9] 
covb7_b4 <- vcov(table5rb)[14,9]
covb7_b5 <- vcov(table5rb)[14,10]
covb2_b0 <- vcov(table5rb)[3,1]
covb3_b0 <- vcov(table5rb)[4,1]
covb6_b0 <- vcov(table5rb)[11,1]
covb3_b2 <- vcov(table5rb)[4,3] 
covb6_b2 <- vcov(table5rb)[11,3]
covb6_b3 <- vcov(table5rb)[11,4]

## stability=2
pdf(file=paste("../Party Politics/graphics/Figure 6.pdf", sep=""), width = 6, height = 6)
z1=0  #supply lower bound for z
z2=5   #supply upper bound for z
ci <- 1.684 #90%
nelections <- 2
z <- seq(z1,6,.1)
fz <- c(z,z)
y1 <- (b1+b4*z+b5*nelections+b7*z*nelections)+ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
y2 <- (b1+b4*z+b5*nelections+b7*z*nelections)-ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
fy <- c(y1,y2)
fline <- (b1+b4*z+b5*nelections+b7*z*nelections)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,
     xlab='Moderator: District Magnitude (Log)',
     ylab='Simple Slope: Ethnic Heterogeneity',
     main='Figure 6. Marginal Effect of Ethnic Heterogeneity on ENPV, \nby DM when Stability=2',
     xlim=c(0,5))
lines(z,fline)
f0 <- array(0,c(length(z)))
lines(z,f0,col=8)
zerocross <- round(uniroot.all(approxfun(z,y1), interval = range(c(0,10))),digits=3) ## check for zero-crossings 
legend("bottomright", cex=.5,  bty = "n", bg="white", legend = "Table 5")
dev.off()

###################################################
# FIGURE 7: NINE (9) ELECTIONS
# 
#
###################################################

table5rb <- lm_robust(enpv_bn ~ ethnic + logdm +  
                        dep_elec_acc +
                        (logdm*ethnic) +
                        (logdm*dep_elec_acc) +
                        (dep_elec_acc *ethnic) +
                        (dep_elec_acc*ethnic) +
                        (logdm*dep_elec_acc*ethnic) +
                        upper + 
                        (upper*ethnic) + 
                        proximity + 
                        enc_bn + 
                        (proximity*enc_bn) + 
                        time,
                      data=bunkernegretto2022,
                      clusters=country_name,
                      se_type="stata")
summary(table5rb)

# get betas
b0  <- table5rb$coefficients[1]   
b1  <- table5rb$coefficients[2]   
b2  <- table5rb$coefficients[3]   
b3  <- table5rb$coefficients[4]   
b4  <- table5rb$coefficients[9]  
b5  <- table5rb$coefficients[10]  
b6  <- table5rb$coefficients[11] 
b7  <- table5rb$coefficients[14]  
# get variance
v0 <- vcov(table5rb)[1,1]
v1 <- vcov(table5rb)[2,2]
v2 <- vcov(table5rb)[3,3]
v3 <- vcov(table5rb)[4,4]
v4 <- vcov(table5rb)[9,9] 
v5 <- vcov(table5rb)[10,10]
v6 <- vcov(table5rb)[11,11]
v7  <- vcov(table5rb)[14,14]
## covariance
covb4_b1 <- vcov(table5rb)[9,2]
covb5_b1 <- vcov(table5rb)[10,2]
covb7_b1 <- vcov(table5rb)[14,2]
covb5_b4 <- vcov(table5rb)[10,9] 
covb7_b4 <- vcov(table5rb)[14,9]
covb7_b5 <- vcov(table5rb)[14,10]
covb2_b0 <- vcov(table5rb)[3,1]
covb3_b0 <- vcov(table5rb)[4,1]
covb6_b0 <- vcov(table5rb)[11,1]
covb3_b2 <- vcov(table5rb)[4,3] 
covb6_b2 <- vcov(table5rb)[11,3]
covb6_b3 <- vcov(table5rb)[11,4]

## stability=9
pdf(file=paste("../Party Politics/graphics/Figure 7.pdf", sep=""), width = 6, height = 6)
z1=0  #supply lower bound for z
z2=5   #supply upper bound for z
ci <- 1.684 #90%
nelections <- 9
z <- seq(z1,6,.1)
fz <- c(z,z)
y1 <- (b1+b4*z+b5*nelections+b7*z*nelections)+ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
y2 <- (b1+b4*z+b5*nelections+b7*z*nelections)-ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
fy <- c(y1,y2)
fline <- (b1+b4*z+b5*nelections+b7*z*nelections)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,
     xlab='Moderator: District Magnitude (Log)',
     ylab='Simple Slope: Ethnic Heterogeneity',
     main='Figure 7. Marginal Effect of Ethnic Heterogeneity on ENPV, \nby DM when Stability=9',
     xlim=c(0,5))
lines(z,fline)
f0 <- array(0,c(length(z)))
lines(z,f0,col=8)
zerocross <- round(uniroot.all(approxfun(z,y2), interval = range(c(0,10))),digits=3) ## check for zero-crossings 
abline(v=zerocross,col=4,lty=2)
legend("topright", cex=.5,  bty = "n", bg="white", legend = paste0("Significant above ", zerocross))
legend("bottomright", cex=.5,  bty = "n", bg="white", legend = "Table 5")
dev.off()

###################################################
# FIGURE A2: ALTERNATIVE PLUG-IN ELECTIONS
# 
#
###################################################

# ## number=variable (plug in): 
# nelections <- 12 ##        <------------------- plug-in number of elections
# pdf(file=paste("../Party Politics/graphics/Figure A2 (",nelections," elections).pdf", sep=""), width = 6, height = 6)
# z1=0  #supply lower bound for z
# z2=5   #supply upper bound for z
# ci <- 1.684 #90%
# z <- seq(z1,6,.1)
# fz <- c(z,z)
# y1 <- (b1+b4*z+b5*nelections+b7*z*nelections)+ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
# y2 <- (b1+b4*z+b5*nelections+b7*z*nelections)-ci*sqrt(v1+(z^2)*v4+((nelections)^2)*v5+(z^2)*((nelections)^2)*v7+2*z*covb4_b1+2*nelections*covb5_b1+2*z*nelections*covb7_b1+2*z*nelections*covb5_b4+2*nelections*(z^2)*covb7_b4+2*((nelections)^2)*z*covb7_b5)
# fy <- c(y1,y2)
# fline <- (b1+b4*z+b5*nelections+b7*z*nelections)
# plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,
#      xlab='Moderator: District Magnitude (Log)',
#      ylab='Simple Slope: Ethnic Heterogeneity',
#      main=paste0("Stability=",nelections),
#      xlim=c(0,5))
# lines(z,fline)
# f0 <- array(0,c(length(z)))
# lines(z,f0,col=8)
# legend("bottomright", cex=.5,  bty = "n", bg="white", legend = "Table 3")
# dev.off()

###################################################
# TABLE A1/MODEL 1
#
# 
###################################################

## other interaction
bunkernegretto2022$sp <- bunkernegretto2022$dep_dm_1tier*bunkernegretto2022$dep_as_1tier
d_a1m1 <- bunkernegretto2022
table_a1m1 <- lm_robust(enpv_bn ~ log(sp) + 
                          upper + pres_formula + proximity + 
                          (pres_formula*proximity) + 
                          time, data=d_a1m1,
                        clusters=country_name,
                        se_type="stata")
summary(table_a1m1)
nobs(table_a1m1)

###################################################
# TABLE A1/MODEL 2
#
# 
###################################################

bunkernegretto2022$sp <- bunkernegretto2022$dep_dm_1tier*bunkernegretto2022$dep_as_1tier
d_a1m2 <- subset(bunkernegretto2022, number < 6)
table_a1m2 <- lm_robust(enpv_bn ~ log(sp) + 
                          upper + pres_formula + proximity + 
                          (pres_formula*proximity) + 
                          time, data=d_a1m2,
                        clusters=country_name,
                        se_type="stata")
summary(table_a1m2)

###################################################
# TABLE A1/MODEL 3
#
# 
###################################################

d_a1m3 <- subset(bunkernegretto2022, number > 5)
table_a1m3 <- lm_robust(enpv_bn ~ log(sp) + 
                          upper + pres_formula + proximity + 
                          (pres_formula*proximity) + 
                          time, data=d_a1m3,
                        clusters=country_name,
                        se_type="stata")
summary(table_a1m3)

###################################################
# TABLE A2/MODEL 1
#
# 
###################################################

bunkernegretto2022$sp <- bunkernegretto2022$dep_dm*bunkernegretto2022$dep_as
bunkernegretto2022$spm <- (((bunkernegretto2022$sp)^(1/4)) + 1)^(2/3)
d_a2m1 <- subset(bunkernegretto2022, upper==0)
table_a2m1 <- lm_robust(enpv_bn ~ spm, 
                        data=d_a2m1,
                        clusters=country_name,
                        se_type="stata")
summary(table_a2m1)
nobs(table_a2m1)

###################################################
# TABLE A2/MODEL 2
#
# 
###################################################

## 
bunkernegretto2022$spm <- ((bunkernegretto2022$sp)^(1/4) + 1)^(2/3)
d_a2m2 <- subset(bunkernegretto2022, upper==0 & number < 6)
table_a1m2 <- lm_robust(enpv_bn ~ spm, 
                        data=d_a2m2,
                        clusters=country_name,
                        se_type="stata")
summary(table_a1m2)
nobs(table_a1m2)

###################################################
# TABLE A2/MODEL 3
#
# 
###################################################

##
bunkernegretto2022$spm <- ((bunkernegretto2022$sp)^(1/4) + 1)^(2/3)
d_a2m3 <- subset(bunkernegretto2022, upper==0 & number > 5)
table_a1m3 <- lm_robust(enpv_bn ~ spm, 
                        data=d_a2m3,
                        clusters=country_name,
                        se_type="stata")
summary(table_a1m3)
nobs(table_a1m3)

###################################################
# TABLE A2/MODEL 4
#
# 
###################################################

## Generate SMP2 (eq. 8: Vote Winning Parties: complex)
j <- 1.7 # empiriclly derrived kb
#j <- 2.5 # s+t number
bunkernegretto2022$msb   <-  bunkernegretto2022$dep_as_1tier * bunkernegretto2022$dep_dm_1tier
bunkernegretto2022$t     <-  bunkernegretto2022$upper
bunkernegretto2022$spm2  <-  ((j)^(bunkernegretto2022$t)*((bunkernegretto2022$msb)^(1/4)) + 1)^(2/3)

d_a2m4 <- subset(bunkernegretto2022)
table_a2m4 <- lm_robust(enpv_bn ~ spm2, 
                        data=d_a2m4,
                        clusters=country_name,
                        se_type="stata")
summary(table_a2m4)
nobs(table_a2m4)

###################################################
# TABLE A2/MODEL 5
#
# 
###################################################

## Generate SMP2 (eq. 8: Vote Winning Parties: complex)
j  <- 1.7 # empirically derived kb
#j <- 2.5 # s+t number
bunkernegretto2022$msb    <- bunkernegretto2022$dep_as_1tier * bunkernegretto2022$dep_dm_1tier
bunkernegretto2022$t      <-  bunkernegretto2022$upper
bunkernegretto2022$spm2  <-  ((j)^(bunkernegretto2022$t)*((bunkernegretto2022$msb)^(1/4)) + 1)^(2/3)

d_a2m5 <- subset(bunkernegretto2022, number < 6)
table_a2m5 <- lm_robust(enpv_bn ~ spm2, 
                        data=d_a2m5,
                        clusters=country_name,
                        se_type="stata")
summary(table_a2m5)
nobs(table_a2m5)

###################################################
# TABLE A2/MODEL 6
#
# 
###################################################

## Generate SMP2 (eq. 8: Vote Winning Parties: complex)
j  <- 1.7 # empirically derived kb
#j <- 2.5 # s+t 
bunkernegretto2022$msb    <- bunkernegretto2022$dep_as_1tier * bunkernegretto2022$dep_dm_1tier
bunkernegretto2022$t      <-  bunkernegretto2022$upper
bunkernegretto2022$spm2  <-  ((j)^(bunkernegretto2022$t)*((bunkernegretto2022$msb)^(1/4)) + 1)^(2/3)

d_a2m6 <- subset(bunkernegretto2022, number > 5)
table_a2m6 <- lm_robust(enpv_bn ~ spm2, 
                        data=d_a2m6,
                        clusters=country_name,
                        se_type="stata")
summary(table_a2m6)
nobs(table_a2m6)
nobs(table_a2m6)

###################################################
# TABLE A3/MODEL 1
#
# 
###################################################

bunkernegretto2022$sp <- bunkernegretto2022$dep_dm_1tier*bunkernegretto2022$dep_as_1tier
d_a3m1 <- subset(bunkernegretto2022, dep_elec_acc>2)
table_a3m1 <- lm_robust(enpv_bn ~ log(sp) + 
                          upper + pres_formula + proximity + 
                          (pres_formula*proximity) + 
                          time, data=d_a3m1,
                        clusters=country_name,
                        se_type="stata")
summary(table_a3m1)
nobs(table_a3m1)

###################################################
# TABLE A3/MODEL 2
#
# 
###################################################

d_a3m2 <- subset(bunkernegretto2022, dep_elec_acc>2 & number < 6)
table_a3m2 <- lm_robust(enpv_bn ~ log(sp) + 
                          upper + pres_formula + proximity + 
                          (pres_formula*proximity) + 
                          time, data=d_a3m2,
                        clusters=country_name,
                        se_type="stata")
summary(table_a3m2)
nobs(table_a3m2)

###################################################
# TABLE A3/MODEL 3
#
# 
###################################################

d_a3m3 <- subset(bunkernegretto2022, dep_elec_acc>2 & number > 5)
table_a3m3 <- lm_robust(enpv_bn ~ log(sp) + 
                          upper + pres_formula + proximity + 
                          (pres_formula*proximity) + 
                          time, data=d_a3m3,
                        clusters=country_name,
                        se_type="stata")
summary(table_a3m3)
nobs(table_a3m3)

###################################################
# TABLE A4/MODEL 1
#
# 
###################################################

bunkernegretto2022$sp <- bunkernegretto2022$dep_dm_1tier*bunkernegretto2022$dep_as_1tier
bunkernegretto2022$spm <- ((bunkernegretto2022$sp)^(1/4) + 1)^(2/3)
d_a4m1 <- subset(bunkernegretto2022, upper==0 & dep_elec_acc>2)
table_a4m1 <- lm_robust(enpv_bn ~ spm, 
                        data=d_a4m1,
                        clusters=country_name,
                        se_type="stata")
summary(table_a4m1)

###################################################
# TABLE A4/MODEL 2
#
# 
###################################################

bunkernegretto2022$spm <- ((bunkernegretto2022$sp)^(1/4) + 1)^(2/3)
d_a4m2 <- subset(bunkernegretto2022, upper==0 & number < 6  & dep_elec_acc>2)
table_a4m2 <- lm_robust(enpv_bn ~ spm, 
                        data=d_a4m2,
                        clusters=country_name,
                        se_type="stata")
summary(table_a4m2)

###################################################
# TABLE A4/MODEL 3
#
# 
###################################################

bunkernegretto2022$spm <- ((bunkernegretto2022$sp)^(1/4) + 1)^(2/3)
d_a4m3 <- subset(bunkernegretto2022, upper==0 & number > 5 & dep_elec_acc>2)
table_a4m3 <- lm_robust(enpv_bn ~ spm, 
                        data=d_a4m3,
                        clusters=country_name,
                        se_type="stata")
summary(table_a4m3)

###################################################
# TABLE A4/MODEL 4
#
# 
###################################################

## Generate SMP2 (eq. 8: Vote Winning Parties: complex)
j  <- 1.7 # empirically derived kb
#j <- 2.5 # s+t number
bunkernegretto2022$msb    <- bunkernegretto2022$dep_as_1tier * bunkernegretto2022$dep_dm_1tier
bunkernegretto2022$t      <-  bunkernegretto2022$upper
bunkernegretto2022$spm2   <-  (((j)^(bunkernegretto2022$t)) * ((bunkernegretto2022$msb)^(1/4)) + 1)^(2/3)

d_a4m4 <- subset(bunkernegretto2022, dep_elec_acc>2)
table_a4m4 <- lm_robust(enpv_bn ~ spm2, 
                        data=d_a4m4,
                        clusters=country_name,
                        se_type="stata")
summary(table_a4m4)
nobs(table_a4m4)

###################################################
# TABLE A4/MODEL 5
#
# 
###################################################

## Generate SMP2 (eq. 8: Vote Winning Parties: complex)
j  <- 1.7 # empirically derived kb
#j <- 2.5 # s+t 
bunkernegretto2022$msb    <- bunkernegretto2022$dep_as_1tier * bunkernegretto2022$dep_dm_1tier
bunkernegretto2022$t2      <-  (bunkernegretto2022$dep_as-bunkernegretto2022$dep_as_1tier)/bunkernegretto2022$dep_as
bunkernegretto2022$t      <-  bunkernegretto2022$upper
bunkernegretto2022$spm2   <-  (((j)^(bunkernegretto2022$t)) * ((bunkernegretto2022$msb)^(1/4)) + 1)^(2/3)

d_a4m5 <- subset(bunkernegretto2022, number < 6 & dep_elec_acc>2)
table_a4m5 <- lm_robust(enpv_bn ~ spm2, 
                        data=d_a4m5,
                        clusters=country_name,
                        se_type="stata")
summary(table_a4m5)
nobs(table_a4m5)

###################################################
# TABLE A2/MODEL 6
#
# 
###################################################

## Generate SMP2 (eq. 8: Vote Winning Parties: complex)
j  <- 1.7 # empirically derived kb
#j <- 2.5 # s+t
bunkernegretto2022$msb    <- bunkernegretto2022$dep_as_1tier * bunkernegretto2022$dep_dm_1tier
bunkernegretto2022$t      <-  bunkernegretto2022$upper
bunkernegretto2022$spm2   <-  (((j)^(bunkernegretto2022$t)) * ((bunkernegretto2022$msb)^(1/4)) + 1)^(2/3)

d_a4m6 <- subset(bunkernegretto2022, number > 5 & dep_elec_acc>2)
table_a4m6 <- lm_robust(enpv_bn ~ spm2, 
                        data=d_a4m6,
                        clusters=country_name,
                        se_type="stata")
summary(table_a4m6)
nobs(table_a4m6)
